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ABSTRACT 


A  model  of  a  microwave-induced  plasma  propulsion  system  has  been 
developed  in  one  dimension  for  a  transverse  electric  mode  (TEio)  of  operation  in  a 
rectangular  waveguide.  Available  experimental  data  are  compared  to  the 
computational  results  for  the  case  of  a  planar  propagating  plasma  wave  and,  using 
a  TEio  mode-shape  approximation,  for  a  wave  propagating  in  a  waveguide. 
Temperature  profiles,  plasma  propagation  velocities,  velocity  profiles,  and  absorbed 
power  histories  are  obtained  for  flow  of  helium  from  .5  to  1  atmosphere  pressure 
and  500  to  3000  watts  input  power  at  a  frequency  of  2.45  GHz.  The  computational 
results  show  the  observed  jumping  of  the  plasma  towards  the  microwave  source. 
Peak  plasma  temperatures  range  from  8000  to  9000  K  over  the  input  power  range. 
For  an  input  power  of  1081.7  W  the  calculated  percentage  of  power  absorbed  is 
approximately  70  percent  for  the  planar  case  and  40  percent  for  the  waveguide  case. 
Comparisons  with  experimental  data  indicate  other  mechanisms  (not  involving 
transient  processses),  most  likely  associated  with  the  nonequilibriiun  behavior  of 
the  plasma,  are  responsible  for  the  disagreement  between  the  model  results  and  the 
observed  plasma  propagation  velocities. 
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1.  INTRODUCTION 


1.1  OVERVIEW  -  MICROWAVE  PROPULSION 

Some  of  the  advanced  propulsion  systems  currently  being  considered  for 
future  space  missions  are  the  resistojet,  arcjet,  ion  engines,  and  laser  and 
microwave  electrothermal  propulsion  systems.  Studies  conducted  of  advanced 
propulsion  concepts  for  both  NASA  and  Air  Force  missions,  such  as  low-earth  orbit 
transfers  and  long  duration  missions,  have  yielded  certain  optimum  parameters  of 
operation.  With  current  electrical  power  generation  capabilities,  propulsion 
systems  with  specific  impulses  in  the  range  of  1000  to  2000  seconds  would  optimi? 
the  payload  ratio  of  these  spacecraft  (Micci,  1984).  Of  the  advanced  propulsion 
concepts  mentioned,  microwave  electrothermal  propulsion  systems  best  suit  this 
range  of  operation. 

Resistojets  employ  electrothermal  heating  through  electrical  resistance  in  a 
solid  placed  in  the  flow  of  the  propellant.  These  devices  are,  therefore,  limited  by 
the  materials  involved  and  currently  achieve  a  maximvun  specific  impulse  of  below 
385  seconds  even  when  using  hydrogen  as  the  propellant  (Mueller  and  Micci,  1989). 
Arcjets  also  use  electrothermal  heating,  but  do  so  through  use  of  an  electrical  arc 
passing  through  the  fluid.  These  systems  have  a  specific  impidse  of  approximately 
700  seconds  for  ammonia  propellants  at  high  input  power  levels  (Mueller  eind  Micci, 
1990).  The  electrodes  are  subject  to  erosion,  especially  at  higher  input  power  levels, 
and  thereby  limit  the  lifetime  of  such  systems.  In  the  study  of  cathode  erosion  by 
Deininger,  Chopra,  and  Goodfellow,  the  lifetime  of  a  30  kW-class  arcjet  was  only 
573  hours  (Deininger,  et  al.,  1989)  as  compared  to  the  lifetime  estimate  of  10,000 
hours  for  such  mission  profiles  as  orbit  transfer  and  orbit  control  for  the  space 
station  (Mueller  and  Micci,  1990).  A  propulsion  system  using  arcjets  wotJd  thus 
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reqmrc  backup  thrusters  adding  to  the  overall  mass  of  the  spacecraft.  Propulsion 
devices  relying  on  electrostatic  and  electromagnetic  means  have  optimum  operation 
at  much  higher  specific  impulses  and  have  thrust  levels  on  the  order  of  only  a  few 
Newtons. 

Most  laser  thermal  propulsion  concepts  use  the  laser  to  optically  heat  the 
propellant  forming  a  plasma,  or  hot  gas  region.  The  heated  propellant  is  then 
accelerated  by  use  of  a  standard  supersonic  nozzle  configuration.  Laser  thermal 
systems  could  have  the  laser  energy  produced  elsewhere,  such  as  on  the  Earth,  and 
beamed  to  the  spacecraft,  making  for  a  very  low  specific  mass.  Problems  arise, 
however,  in  focusing,  absorption,  refraction,  and  divergence  of  the  beam  through  the 
atmosphere.  Also,  laser  systems  can  suffer  from  low  efficiencies  due  to  radiation 
losses  from  the  plasma  and  low  energy  conversion  processes  (Jeng  and  Keefer, 
1987).  Microwave  thermal  propulsion  devices  are  also  an  electrodeless,  plasma 
producing  device  and  have  better  efficiencies  than  those  of  laser  systems.  Gas 
absorptivities  are  higher  at  microwave  wavelengths  than  at  optical  wavelengths  for 
temperatures  below  those  yielding  significant  ionization.  This  enables  operation  at 
lower  power  levels  and  thus  lower  temperatures,  thereby  eliminating  the  problems 
of  radiative  and  ionization  losses.  Also,  microwave  conversion  efficiencies  are 
higher  and  the  microwave  generating  devices  are  less  complex  (Venkateswaran,  et 
al.,  1990). 

Microwaves  are  generally  considered  to  have  frequencies  ranging  from  10®  to 
1012  Hz,  corresponding  to  characteristic  wavelengths  of  .3  to  .0003  meters.  Since 
these  wavelengths  can  be  on  the  same  order  as  propulsion  chamber  dimensions, 
resonance  characteristics  can  be  obtained  by  proper  dimensioning  of  an  enclosed 
waveguide  shape.  This  is  the  general  reasoning  behind  the  use  of  resonant  cavity 
waveguides.  Microwaves  are  governed  by  Maxwell's  equations  and  as  such  have  an 
infinite  number  of  different  possible  solutions,  each  of  which  corresponds  to  a  mode 
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of  propagation  in  a  waveguide.  These  modes  are  divided  into  three  classifications 
according  to  the  fields  present  in  the  waveguide  when  that  mode  is  propagating. 
These  classifications  are  the  TE,  TM,  and  TEM  modes.  In  the  TE,  or  transverse 
electric,  mode  there  is  no  longitudinal  component  of  the  electric  field,  but  only 
transverse  components  of  the  electric  field.  In  the  TM,  or  transverse  magnetic, 
mode  there  is  no  longitudinal  component  of  the  magnetic  field.  A  wave  with  only 
transverse  components  of  both  the  electric  and  magnetic  fields  is  called  transverse 
electric  and  magnetic,  or  TEM.  The  TE  and  TM  modes  are  further  divided  by  the 
values  of  two  integers,  m  and  n,  used  in  the  expressions  for  the  components  of  the 
field.  These  integers  appear  as  subscripts  of  TE  or  TM.  Different  modes  of 
operation  have  different  frequencies  below  which  there  is  no  propagation  of  the 
microwave.  This  is  called  the  cut-off  condition  (Fuller,  1969).  The  lowest,  or 
dominant,  mode  of  operation  is  usually  used  in  microwave  systems  through  use  of  a 
low  input  frequency.  This  avoids  multiple  modes  of  operation  of  the  microwave 
fields  through  the  use  of  the  cut-off*  condition. 

Several  concepts  for  microwave  propulsion  systems  have  been  considered  in 
the  past.  Micci  (1984)  groups  the  various  systems  into  five  different  categories.  In 
the  first,  the  propellant  is  passed  through  a  constraining  tube  located  within  a 
waveguide  or  resonant  cavity.  This  concept  suffers  from  a  high  temperature  located 
close  to  the  walls  of  the  constraining  tube  which  leads  to  thermal  breakdown.  Also, 
it  was  found  that  heat  loss  by  conduction  to  the  walls  of  the  tube  was  greater  than 
that  of  the  heat  convected  to  the  propellant,  making  such  a  system  undesirable 
(Micci,  1984).  A  second  setup  is  referred  to  as  the  "coaxial  microwave  plasmatron". 
As  the  name  suggests,  the  microwave  energy  is  introduced  through  the  use  of  a 
coaxial  line  instead  of  a  waveguide.  The  propellant  is  introduced  upstream  of  the 
termination  of  the  center  conductor  of  the  coaxial  line  and  the  plasma  is  formed  at 
the  end  of  the  coaxial  line.  Difficulties  arise  from  this  concept  due  to  the  plasma 
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region  being  next  to  the  center  conductor  (Batenin,  et  al.,  1976a;  Batenin,  et  al., 
1975;  Batenin,  et  al.,  1976c;  Miyake,  et  al.,  1974). 

The  third  and  fourth  concepts  involve  the  use  of  a  resonant  cavity 
configuration  in  a  cylindrical  waveguide.  In  the  filamentary  mode,  the  microwave 
energy  is  input  at  the  TMqi  mode  of  operation.  In  this  mode  shape,  the  peak  of  the 
electric  field,  and  thus  the  location  of  the  plasma,  lies  along  the  centerline  of  the 
cylinder.  Experimental  studies  of  this  system  were  done  in  the  Soviet  Union  by 
Kapitza  (1970, 1972)  and  later  by  Tishchenko  and  Zatsepin  (1975)  and  in  Japan  by 
Hasegawa,  et  al.  (1975).  Theoretical  modeling  was  done  by  Mieroivich  (1972a, 
1972b,  1973)  and  by  Asinovskii  and  Batenin  (1973).  For  the  fourth  concept,  the 
TEoi  mode  of  operation  is  used  instead  of  the  TMoi.  In  this  mode  shape,  the  electric 
field  lines  are  circular  and  parallel  to  the  longitudinal  axis  of  the  resonant  cavity 
cylinder.  The  maximiun  of  the  electric  field  occurs  in  a  ring  lying  away  from  the 
walls  of  the  cavity.  The  plasma  will  then  form  in  a  toroidal  shape  along  the  ring  of 
maximum  electric  field  strength.  Further  theoretical  and  experimental 
investigation  of  resonant  cavities  for  use  in  microwave  propulsion  systems  have 
been  conducted  by  Whitehair  and  Asmussen  (1984)  and  Whitehair,  et  al.  (1987)  and 
later  by  Mueller  and  Micci  (1990)  and  Venkateswaran,  et  al.  (1990). 

The  fifth  concept  involves  the  absorption  of  the  microwave  energy  by  a 
plasma  region  contained  within  a  waveguide  serving  as  a  propulsion  chamber.  In 
this  scheme  the  plasma  propagates  toward  the  source  of  the  microwave  energy  and 
must  be  stabilized  by  some  means.  The  plasma  front  then  acts  similarly  to  a 
combustion  wave  front,  heating  the  propellant  which  is  accelerated  out  a  nozzle. 
Experimental  and  theoretical  studies  relevant  to  this  concept  have  been  conducted 
by  Beust  and  Ford  (1961),  Raizer  (1972),  Batenin,  et  al.  (1973, 1976, 1977),  Knecht 
and  Micci  (1988),  Mueller  and  Micci  (1988, 1989, 1990).  The  research  described  in 
this  thesis  is  of  an  alternate  approach  to  the  computational  modeling  of  the  plasma 
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wavefront  concept  and  shows  the  results  and  comparisons  of  the  model  to  available 
experimental  data. 

1.2  OVERVIEW  -  PLASMA  DEPOSITION 

Microwave  induced  plasmas  are  also  ideal  for  coating  applications  such  as 
chemical  vapor  deposition  (CVD),  thin-film  deposition  and  etching  (Asmussen,  1989; 
Hsu,  et  al.,  1989).  Due  to  their  electrodeless  nature,  the  plasma  cannot  be 
contaminated  by  impurities  contained  in  the  electrode  material.  Thus  a  higher 
degree  of  flexibility  in  the  composition  of  the  plasma  is  possible  with  respect  to  the 
radiation  properties.  This  is  of  great  importance  in  CVD  processes.  Microwaves 
allow  for  better  control  of  the  temperature  distribution  through  variation  of  the 
input  frequency.  Also,  manufacturing  and  lifetime  aspects  are  improved  through 
the  use  of  microwave  systems  (Offermanns,  1989).  Further  understanding  of  the 
d3mamics  of  the  interactions  of  the  gas  or  gases  with  the  plasma  would  be  beneficial 
in  the  optimization  of  various  coating  qualities.  The  following  analysis  and  results, 
however,  focus  on  the  propulsion  aspect  of  microwave  induced  plasmas,  although 
the  model  could  be  similarly  applied  to  model  plasmas  used  in  chemical  vapor 
deposition  applications. 

1.3  MODELING  APPROACH 

A  microwave  discharge  in  a  waveguide  is  observed  to  propagate  towards  the 
source  of  the  power  in  configurations  with  no  added  stabilization  techniques.  There 
are  four  main  propagation  mechanisms  that  may  be  involved.  Three  mechanisms 
described  in  the  paper  by  Bethke  and  Reuss  (1969)  for  low  pressures  are  gas 
photoionization  due  to  resonant  radiation,  electron  ionization  due  to  collisions  with 
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energetic  electrons  and  electron  difEusion.  Myshenkov  and  Raizer  (1971)  combined 
these  first  two  mechanisms  as  resonant  radiation.  For  pressures  closer  to 
atmospheric,  the  dominant  mechanism  is  heat  conduction  (Batenin,  et  al.,  1973; 
Raizer,  1972).  Batenin,  et  al.  (1973)  also  foimd  resonant  radiation  important  for 
heavy  gases  at  atmospheric  pressures.  In  each  case,  the  particular  mechanism 
causes  microwave  heating  upstream  of  the  existing  plasma  region  causing  the 
plasma  to  move  towards  the  source  of  the  power  and  extinguishment  of  the  plasma 
downstream.  The  extinguishment  develops  due  to  attenuation  of  the  microwaves 
through  the  plasma  region,  and  thus,  the  downstream  region  cannot  maintain  the 
plasma. 

Theoretical  models  of  microwave  systems  have  been  developed  using  several 
different  approaches.  Analytic  solutions  of  the  electromagnetic  system  have  been 
developed  by  Myshenkov  and  Raizer  (1972)  and  Raizer  (1972).  In  each,  the 
equations  are  solved  in  the  frequency  domain,  and  therefore,  assume  a  frequency 
dependence  of  the  electric  held.  Several  simplifying  assiunptions  are  used  such  that 
the  particular  problem  is  reduced  to  an  eigenvalue  problem  for  the  propagation 
velocity. 

In  computational  models  of  electromagnetic  field  propagation,  there  are  two 
principle  approaches  to  the  problem.  In  the  first  approach,  the  electromagnetic 
system  is  solved  in  the  time  domain.  This  puts  severe  restrictions  on  the  maximum 
allowable  time  step.  Yee  (1966)  developed  an  algorithm  for  solutions  of  Maxwell's 
equations.  This  model  was  later  improved  by  Mur  (1981)  who  modified  the 
boundary  conditions,  and  by  Goorjian  (1990).  The  work  of  Venkateswaran,  Merkle, 
and  Micci  (1990)  is  based  on  the  algorithm  developed  by  Yee. 

In  the  second  modeling  approach,  the  equations  are  solved  in  the  frequency 
domain,  and  thus,  the  electric  field  has  an  assumed  frequency  dependence.  This 
approach  eliminates  the  problems  of  the  time  increment.  Rhodes  and  Keefer  (1989) 
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developed  an  axisymmetric,  two  dimensional  model  for  a  radio  frequency  gas  heater 
which  can  be  related  to  a  microwave  system.  They  used  the  frequency  domain  and 
solved  the  system  by  use  of  a  complex,  electric  field  vector  potential  in  the  SIMPLE 
algorithm  derived  finm  earlier  work.  In  the  works  of  Sforza  (1969)  and  Ofiermanns 
(1989),  models  of  the  microwave  system  were  developed  using  series  expansion 
approximations. 

The  objectives  of  this  research  include  formulation  of  an  electromagnetic 
equation  system  coupled  to  a  full  time  dependent  fluid  solution  to  better  model  the 
interactions  observed  in  microwave  electrothermal  propulsion  devices.  Among  these 
improvements  is  the  need  of  a  model  that  will  show  the  observed  phenomenon  of 
"hopping"  of  the  plasma  within  the  waveguide.  The  anomalous  values  for  power 
absorption  between  the  model  of  Mueller  and  Micd  (1989)  with  experimental  data 
(Batenin,  et  al.,  1976b;  Mueller  and  Micci,  1989;  Raizer,  1972)  must  be  resolved. 
Also,  it  is  desired  to  obtain  the  optimal  input  parameters  that  will  result  in  peak 
input  power  absorption  by  the  plasma. 
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2.  THEORETICAL  ANALYSIS 


2.1  FORMULATION  OF  THE  GOVERNING  EQUATIONS 


The  basic  flow  geometry  in  a  waveguide  of  either  a  rectangular  or  circular 
cross  section  is  shown  schematically  in  Figures  la,  lb,  and  2.  The  plasma  front  is 
assumed  to  be  initiated  within  the  waveguide  and  will  propagate  toward  the 
microwave  source  unless  some  means  of  stabilization  is  provided.  Thus,  we  consider 
a  non-steady  formulation  of  the  fluid  dynamic  equations.  Since  the  considered 
operating  pressure  is  approximately  atmospheric  and  flow  or  propagation  speeds 
are  low  •  subsonic,  the  plasma  is  assumed  to  be  in  local  thermodynamic  equilibrium. 
Additionally,  due  to  absorption  limitations,  the  plasma  ionization  fraction  is 
typically  less  than  0.01 ,  and  radiative  transfer  is  assumed  negligible  based  on  the 
results  of  Mueller  and  Micd  (1989). 


Figure  la.  Coordinate  convention  for  a  rectangular  wavegfuide. 


8 


z 

Figure  lb.  Coordinate  convention  for  a  cylindrical  waveguide. 


Figure  2.  Schematic  of  plasma  propagation  in  a  rectangular  waveguide. 


The  flow  is  therefore  governed  by  the  Navier-Stokes  equations  for  a 
compressible  gas  in  equilibrium  (see,  for  example,  Pai,  1962) 

^  +  V-  (pu)  =0 

at  ^  (1) 
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(2) 


+  V-  (puu-t)  =noExH 

a  t 


-»2 


+  V-  (puh  +  ^-  u^  =aE 
d  t 


(3) 


together  with  the  thermal  and  caloric  equations  of  state.  The  time-averaged 
Lorentz  force  and  Ohmic  heating  terms  in  Equations  (2)  and  (3)  are  determined  by 
the  Maxwell  equations  given  below. 


dpH 

dt 


-V  X  E 


(4) 


^+aE- VxH  =  0 

at 


(5) 


The  above  equation  uses  Ohm's  Law  to  give  the  current  density  in  terms  of  the 
conductivity  and  the  electric  field. 

For  the  problem  of  interest,  a  single  frequency  field  is  assumed  to  be  output 
by  the  microwave  generator.  The  mediuim  is  assumed  to  be  unmagnetized  such  that 
the  permeability  may  be  replaced  by  the  permeability  of  free  space  and  is  also 
assumed  to  be  free  of  space  charge.  Both  the  inverse  frequency  and  wave  reflection 
time  scales  of  the  electromagnetic  field  are  readily  shown  to  be  much  shorter  than 
any  of  the  fluid  time  scales.  Accordingly,  a  complex  representation  of  the 
electromagnetic  field  in  the  frequency  domain  is  adopted,  such  that 

E(x,t)  =  RcKErCx)  +  iEi(x))e’®*]  (6) 

with  a  corresponding  description  for  the  magnetic  field,  H.  For  the  Transverse 
Electric  (TE)  mode  of  operation,  in  which  the  axial  component  of  the  electric  field  is 
zero  everywhere,  Maxwell's  equations  are  solved  for  the  spatial  amplitude  of  the 
electric  field,  E.  For  the  Transverse  Magnetic  (TM)  mode.  Maxwell's  equations 
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would  be  solved  for  the  spatial  amplitude  of  the  magnetic  field,  H.  The  following 
derivation  is  for  the  TE  mode  of  operation. 

A  multiple  scale  analysis  of  the  equation  was  done  to  determine  the  relevant 
terms.  Two  times,  corresponding  to  the  slow  time  associated  with  the  fluid  time 
scale,  t,  and  fast  time  associated  with  the  frequency  of  the  input  electric  field,  x, 
were  employed.  Thus,  the  value  of  t/x  relevant  to  the  frequencies  used  was  on  the 
order  of  10-®.  Derivatives  with  respect  to  time  of  e  and  a  were  found  to  be  negligible 
in  comparison  to  the  remaining  terms.  Only  zeroth  order  terms  are  retained. 
Splitting  the  resulting  equation  into  real  and  imaginary  parts  yields  coupled 
equations  for  the  real  and  imaginary  spatial  amplitudes,  Ej{  and  Ej.  Since  we  are 
dealing  with  the  TE  mode  of  operation,  however,  there  is  no  axial  (z)  component  of 
the  electric  field.  Thus,  the  equations  may  be  written  as  follows. 

Real:  V^ER(x,y))  +  a)^pe(ER(x,y))  +  (i)pa(Ei(x,y))  *  0  (7) 

Imaginary:  ^  (Ei(x.y))  +  (o^pe (Ei(x,y))  -  (opc  (ER(x.y))  =  0  (8) 

Further  reduction  of  the  vectors  into  scalar  quantities  is  done  in  both  the  one  and 
two  dimensional  formulations  and  shall  be  discussed  in  these  separate  sections. 

2.2  BOUNDARY  CONDITIONS 

The  boundary  conditions  at  the  inlet  and  exit  of  the  waveguide  for  both  the 
one  and  two  dimensional  formulations  make  use  of  the  left  and  right  travelling 
wave  characteristic  equations.  In  the  two  dimensional  formulation,  boundary 
conditions  for  the  wall  and  centerline  are  also  necessary.  These  boundary 
conditions  will  be  discussed  in  the  two  dimensional  section. 


11 


2.2.1  Inlet  Boundary  Conditions 


The  electric  field  at  the  inlet  can  be  expressed  as  the  sum  of  the  known  input 
electric  field,  R,  which  is  a  right  running  wave  and  an  unknown  reflected  electric 
field,  L,  which  is  left  running.  The  characteristic  equations  for  a  left  and  right 
travelling  wave  are  given  below. 

dL  m3L  Q 

at  paz  (9) 

at  p  a z 

The  parameter,  (3,  is  the  phase  constant  (to  be  defined  subsequently)  and  is 
dependent  on  electromagnetic  parameters  and  mode  shape.  The  input  electric  field 
may  be  written  as  follows. 


In  the  two  dimensional  case,  Ro  is  a  function  of  the  radius.  Solving  for  L  in  terms  of 
R  and  (Er  +  iEi)  and  substituting  into  the  left  running  wave  characteristic  equation 
yields 

i  PCEr  +  iEi)  -  =  i2RoPe-*P2 

.  (12) 

At  the  inlet,  the  conductivity  is  assiuned  to  be  zero  so  there  is  no  attenuation  of  the 
wave.  This  assumption  holds  as  long  as  the  inlet  temperature  stays  below 
approximately  3000  K.  The  angular  frequen<y,  ©,  is  assumed  to  be  a  known  input. 
The  amplitude  of  the  input  electric  field,  Ro  (V/m),  is  dependent  on  the  analytic 
solution,  the  averaged  input  power,  and  the  waveguide  dimensions  which  are  all 
related  by  the  power  density  equation: 


(13) 


^input  ^  ®  ^  ^ 

The  term  A  c  refers  to  the  area  of  the  waveguide  cross  section.  Expressions  for  the 
magnetic  field,  H,  in  terms  of  the  electric  field,  E,  are  developed  from  Maxwell's 
equations  and  used  in  the  Equation  (13). 

2.2.2  Exit  Boundary  Conditions 

At  the  exit  a  non-reflective  boxmdary  condition  is  employed.  Here  the  electric 
field  must  satisfy  the  characteristic  equation  for  a  right  running  wave. 

a  (Er  +  iEj)  ^  ^  (Er  iEi)  ^  q 

at  P  az  (14) 

It  should  be  noted,  however,  that  once  the  plasma  has  formed,  there  is  no 
transmittance  of  the  wave  past  the  exit  due  to  attenuation  of  the  wave. 

2.3  COUPLING  TERMS 


For  the  energy  equation,  the  coupling  involves  the  term  E  .  Since  the 
frequency  domain  is  employed,  the  Ohmic  heating  term  must  be  averaged  over  a 
period  of  oscillation.  Returning  to  the  time  dependent  mode  and  averaging  over  a 
period  of  oscillation,  it  can  be  shown  for  both  the  one  and  two  dimensional  cases 
that 

Z2  1  2  -  2\ 

E=11Eb+E,) 

The  coupling  in  the  momentum  equation  is  due  to  E  x  H.  Using  Maxwell's 
equations,  relations  for  the  magnetic  field  in  terms  of  the  electric  field  components 
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are  developed.  The  cross  product  is  averaged  over  a  period  of  oscillation  like  that  of 
the  £2  term.  The  resulting  equation  is  different  for  the  one  and  two  dimensional 
formulations  and  so  will  be  described  in  each  section. 

2.4  ONE  DIMENSION 


A  one  dimensional  analysis  can  be  used  to  approximate  the  solution  in  a 
rectangular  waveguide  for  the  lowest  TE  mode  of  operation,  the  TEio  mode.  For 
this  case,  the  only  non-zero  component  of  the  electric  field  is  the  y  component  (see 
Figure  1).  Thus,  Er  and  Ei  become  scalar  quantities  representing  the  real  and 
imaginary  y-components  of  the  electric  field.  In  matrix  form  the  coupled  equations 
may  be  written  as  follows. 


-?_  +  (oV 

,dz  I 


(opo 


-  copa 


’ 

+  co^pe 

; 


Er 

El  . 


(16) 


2.4.1  Boundary  Conditions 

The  simplest  form  of  a  one  dimensional  propagating  wave  is  that  of  a  plane 
wave.  For  a  propagating  plane  wave,  the  incoming  electric  field  is  given  by 

R  =  Ae^®^’P^)  (17) 

where  A,  the  amplitude  (V/m),  is  a  real  constant.  The  phase  constant,  p,  is  given  by 
the  relation 

(18) 
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Using  Maxwell's  equations,  an  expression  for  the  x-component  of  the  magnetic  field 
may  be  developed  in  terms  of  R. 

H=— 

V^TTe  (19) 

Using  the  power  density  expression  in  Equation  (13)  and  the  dimensions  a  and  b  for 
the  width  and  height  of  the  waveguide,  the  expression  for  A  becomes 

^input  _  1 

a  b  2  Vp/e.  (20) 


The  inlet  boundary  condition  then  employs  these  parameters  in  Equation  (12). 

For  a  TEio  mode  wave  propagating  in  a  rectangular  wavegmde,  R  is  given  by 
the  following  expression  (Fuller,  1969). 


R  =  A  [<o  VJH  ^  cos(^)J  e^(^^ 


This  two  dimensional  form  is  used  in  Equation  (13)  by  taking  the  peak  value  of  the 
cosine  term.  In  this  derivation  it  is  therefore  assumed  that  the  wave  maintains  the 
sinusoidal  shape  across  the  x  cross  section.  The  two  dimensionality  also  effects  the 
waveguide  phase  constant,  P,  and  produces  a  cut-off  condition,  below  which  no  wave 
is  propagated  down  the  waveguide.  In  general,  for  a  TEmn  mode,  P  is  given  by 


where  a  is  the  x-dimension  and  b  is  the  y-dimension  of  the  waveguide  (Fuller,  1969). 
For  the  TEio  mode,  p  thus  becomes 
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Using  the  assumed  sine  wave  shape,  the  second  derivative  of  x  restating  from  the 
Laplacian  in  the  two  dimensional  equations  can  be  replaced  with  a  constant  value. 


dx" 


E  =  -(ffE 


(24) 


Thus,  Equation  (24)  alters  the  original  equation  matrix  given  in  Equation  (16). 
.2 


— -ffflV-fe) 


\3  z 


(opc 


copa  — -  +  coV-(f)‘ 


2 


/J 


Er 

Ei 


0 


(25) 


Similar  to  the  plane  wave  propagating  case,  an  expression  may  be  developed  for  the 
amplitude  of  the  magnetic  field,  H,  in  terms  of  the  known  input  electric  field 
amplitude,  R.  The  value  of  A  may  then  be  evaluated  using  the  expression  for  power 
density  in  Equation  (14). 


Pinput 

ab  2pc 


0)  Vpe  p  si  cos2(J^) 


(oVpFpsil 


n 


2  2J 


(26) 


As  mentioned  in  the  general  discussion,  the  exit  boundary  for  both  the  planar  and 
waveguide  solutions  employs  the  characteristic  equation  for  a  right  travelling  wave 
to  ensure  a  non-refiective  condition. 


2.4.2  Coupling  Terms 

Recall  that  the  coupling  to  the  energy  equation  is  given  by  Equation  (15).  For 
the  one  dimensional  case,  the  magnetic  field  consists  of  an  x-component  and  a  z,  or 
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axial,  component.  Thus,  the  cross  product  used  in  the  power  density  relation 
becomes 


ExH=^(ErHr  +  EiHi) 


(27) 


where  the  terms  Hr  and  Hi  refer  to  the  x-component  of  the  magnetic  field.  Using 
Maxwell's  equations,  these  magnetic  field  components  may  be  expressed  in  terms  of 
the  electric  field  components  as  follows. 


dz 


Hi=^^ 

^z 


(28a,  28b) 


2.5  TWO  DIMENSIONS  •  AXISYMMETRIC  WAVEGUIDE 


To  solve  the  equation  system  in  a  circular  waveguide,  a  two  dimensional 
axisymmetric  analysis  is  posed.  The  lowest  TE  mode  of  operation  which  has  an 
axisyrametric  propagation  is  the  TEoi  mode.  Recall  that  by  definition  there  is  no 
axial  component  of  the  electric  field.  Since  the  solution  is  assumed  to  be 
axisymmetric,  there  can  be  no  angular  dependence  for  either  the  radial  or  angular 
electric  field  components.  Furthermore,  the  radial  component  can  be  shown  to  be 
zero  by  the  following  Maxwell’s  relation. 

V.  (eE)  =  l^^^  =  0 

^  (29) 

Thus,  the  radial  component  of  the  electric  field  must  be  a  constant  and,  due  to  the 
centerline  condition,  the  constant  is  zero  (see  Figure  lb). 

In  the  matrix  equations,  therefore,  Er  and  Ei  are  scalar  quantities  referring 
to  the  angular  component  of  the  electric  field. 
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Real: 


Imaginary: 


iALALJ_ 

^arl  BtI  5^2 


l^\r±U± 


’“arl  arj 


az 


(Er)  +  a)2^£(ER)  +  a)^a(Ei)  =  0 


(Ej)  +  ©2)i£(Ei)  -  cdho(Er)  =  0 


(30) 


(31) 


For  values  of  small  r  the  radial  derivative  may  lead  to  computational  error  due  to 
the  1/r  term.  The  problem  with  small  values  of  r  can  be  avoided  by  assuming  a  non- 
radial  characteristic  of  the  equations,  and  approximating  the  input  electric  field  by 
either  a  S3rmmetrically  shaped  wave  (sine  wave)  or  by  using  this  sine  wave 
approximation  for  small  values  of  r.  The  dependent  term  in  this  approximation  is 
then  y  instead  of  r  and  the  equations  become  similar  to  the  two  dimensional 
equations  for  a  rectangular  waveguide.  For  the  chosen  axisymmetric  propagation 
mode  (TEoi),  a  sine  wave  approximates  the  analytic  solution  very  well.  An 
approximation  factor  similar  to  that  involved  in  the  1-D  solution  must  be  added  to 
maintain  the  proper  propagation  characteristics.  The  exact  character  will  be 
discussed  along  with  the  boundary  conditions. 


2.5.1  Boundary  Conditions 


For  a  circular  waveguide,  the  incoming  TEoi  electric  field  is  given  by 
Equation  (31) 


R  =  A 


.i(ci)t  -Pz) 


(32) 


where  Jq'  is  the  derivative  of  the  Bessel  function  Jo  with  respect  to  r,  and  a  is  the 
waveguide  radius  (Fuller,  1969).  For  this  case,  the  phase  constant  is  given  by  the 
following  equation. 
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p=V 


C0^^£  -|2^j 


(33) 


As  mentioned,  the  Bessel  fimction  can  be  approximated  by  a  sine  wave  of  amplitude 
27.86. 


R  =  A  [o)  Vp£  ^ (27.86  sii^^l  e -  P  z) 


(34) 


Since  the  Laplacian  of  this  function  does  not  lead  to  and  thus,  the  original 
equation  set  is  not  satisfied,  a  correction  term  must  be  added  to  account  for  the 
difference.  The  matrix  equations  then  become  as  follows. 


real: 


Imaginary: 


2  2 

d  d 

— 2*^ — 2 
3  y  3  z 


2  2 

3  3 

2  _  2 


(Er)  +  o)^|J£(Er)  +  a)pa(Ei)  =  0 


3  y  3  2  J 


(Ei)  +  o)^pe(Ei)  -  (iHxa(ER)  =  0 


(35) 


(36) 


The  relation  for  A  from  the  power  density  equation  then  becomes 


input 


7c  a 


:-A_ 

2pc 


(37) 


As  in  the  1-D  rectangular  case,  these  parameters  are  used  in  Equation  (12)  for  the 
inlet  boxindary  condition.  Also,  the  exit  boundary  uses  the  characteristic  equation 
for  a  right  travelling  wave  given  in  Equation  (14)  and  thereby  only  accounts  for 
axial  propagation  out  of  the  waveguide. 

In  the  two  dimensional  case  it  is  necessary  to  pose  boundary  conditions  for 
the  wall  and  the  centerline  of  the  waveguide.  For  the  method  developed,  these 
boundary  conditions  are  very  simple.  Using  the  assumption  that  the  wall  is  a 
perfectly  conducting  surface,  the  electric  field  must  be  zero  at  the  wall.  Also,  the 
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character  of  the  input  electric  field  in  the  axisymmetric,  TEqi  mode  places  the 
restriction  that  the  electric  field  is  zero  at  the  centerline. 


2.5.2  Coupling  Terms 

As  in  the  case  of  the  1-D  formulation,  the  energy  equation  coupling  term  is 
given  by  Equation  (15).  For  the  two  dimensional,  axisjrmmetric  case,  the  magnetic 
field  consists  of  an  angular  component  and  an  axial  component.  Thus,  the  cross 
product  used  in  the  power  density  relation  has  both  an  axial  component  and  a  radial 
component. 

E^,=^(EbH,r+EiH,i) 

E  xH,  =l{EBHjii  +  EiH,]) 


The  terms  HrR  and  Hri  refer  to  the  radial  component  of  the  magnetic  field  while  the 
terms  HzR  and  Hzi  refer  to  the  axial  components.  Using  Maxwell’s  equations,  these 
magnetic  field  components  may  be  expressed  in  terms  of  the  electric  field 
components  as  follows. 


H  -  1 

HrR  - 


dz 


a  El 


li®  ay 
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_  1 
"’li®  az 


**  1 1 


aEi 


4©  dy 


(40a,  40b) 

(41a,  41b) 


2.6  COMPUTATIONAL  MODEL 


The  equations  for  the  one  dimensional  formulation  posed  in  the  previous 
sections  are  solved  implicitly  using  central  differences  for  the  complex  amplitudes. 
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Er  and  Ej,  by  matrix  inversion  over  the  grid  domain.  This  approach  eliminates  the 
need  for  iteration  of  the  final  electromagnetic  solution  at  each  time  interval.  The 
local  values  of  the  electric  field  are  then  used  explicitly  in  the  unsteady, 
compressible  Navier-Stokes  equations  for  equilibrium  fiow.  The  fluid  solution  is 
found  using  the  one-dimensional  version  of  the  time  accurate,  finite  voliime  method 
developed  by  Ridder  and  Beddini  (1989).  The  fluid  used,  helium,  is  assumed  to 
behave  as  a  perfect  gas  in  local  thermodynamic  equilibrium  (LTE).  This  modeling 
system  sets  no  restrictions  on  the  location  of  the  plasma.  There  is  no  Courant- 
Freidrichs-Lewy  (CFL)  ntunber  restriction  in  the  normal  sense.  The  time  step  is 
limited  only  by  trimcation  error  involved  in  resolving  the  scales  of  interest 
(principally  the  plasma  wave  thickness  and  propagation  speeds),  and  thus  the  time 
step  may  utilize  a  CFL  number  of  order  lO^  or  more  for  the  posed  problem.  This  is  a 
great  advantage  over  schemes  solved  in  the  time  domain  which  must  use  time 
increments  corresponding  to  the  frequency  of  microwave  oscillation  which  is  much 
smaller  than  any  of  the  fluid  time  scales.  The  grid  is  restricted  by  the  thickness  of 
the  plasma  wavefront.  The  grid  cells,  therefore,  have  a  maximum  increment  beyond 
which  the  electric  field  is  not  accurately  defined.  For  the  range  of  inputs  used  in 
this  work,  the  time  increment  was  10-3  seconds  and  the  grid  spacing  was  lO'® 
meters,  corresponding  to  a  CFL  number  of  approximately  600.  In  this  way  the 
plasma-fluid  system  is  modeled  by  a  time  accurate  approach  that  allows  for  easy 
variation  of  the  input  parameters  to  study  the  effects. 
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3.  RESULTS.ANP  PISCUSSIPN 

The  conditions  and  parameters  of  the  present  investigation  are  chosen  to 
correspond  with  the  prior  experiments  of  Batenin,  et  al.  (1977)  and  Mueller  and 
Micci  (1989).  Both  studies  used  similar  arrangements,  although  the  specific 
conditions  of  the  latter  experiment  are  for  a  TEio  mode  at  2.45  GHz.  The  waveguide 
cross  section  is  7.214  x  3.404  cm^.  Both  experiments  also  used  helium  as  a  test 
medium,  and  an  inlet,  cold  gas  temperature  of  300  K  is  assumed. 

As  indicated  by  the  classical  approximate  analysis  of  Raizer  (1972),  the 
plasma  propagation  and  absorption  conditions  are  strongly  dependent  on  transport, 
thermodynamic  and  electrodynamic  properties.  A  literature  review  was  conducted 
to  determine  the  data,  or  the  most  appropriate  methods  for  developing  data,  for 
helium  properties  over  the  approximate  pressure  range  of  0.5  to  1.0  atmospheres 
and  temperature  range  of  300  to  10000  K.  Values  of  parameters  such  as  viscosity, 
thermal  conductivity,  permittivity,  and  conductivity  are  discussed  and  graphically 
presented  in  the  Appendix. 

With  one  exception  that  shall  be  discussed  subsequently,  all  calculations  to  be 
presented  involve  artificially  "igniting"  the  plasma  near  the  exit  of  the  waveguide  by 
providing  an  enhanced  electrical  conductivity  corresponding  to  a  temperature  of 
6800  K,  The  plasma  wave  is  then  monitored  as  it  moves  to  the  inlet  of  the 
waveguide  through  essentially  stagnant  gas.  The  length  of  the  waveguide  used  is 
one  meter,  thereby  ensuring  adequate  monitoring  distances  of  several  free  space 
wavelengths  (~0.2  m).  The  niunber  of  axial  nodes  employed  in  the  calculations  was 
nominally  1000,  and  preliminary  sensitivity  studies  were  performed  to  ensure  that 
this  number  was  svdficient  to  accurately  resolve  the  thin  plasma  wave  thickness  (~5 
cm)  obtained  at  higher  input  power  levels  (3  kW).  The  calculations  shown  represent 
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the  solutions  for  a  planar  propagating  wave  and  for  a  TEio  mode  wave  propagating 
in  a  waveguide.  Recall  that  the  cross  sectional  characteristics  in  the  waveguide 
were  assumed  to  remain  the  same  as  in  the  analytic  solution  for  propagation  in  a 
waveguide  containing  no  plasma. 

Figure  3  shows  the  variation  with  time  of  the  plasma  propagation  velocity. 
The  seeming  thickness  of  the  lines  is  due  to  the  interpolation  scheme  used  to  find 
the  plasma  wavefront.  The  sinusoidal  shape  shows  the  effects  of  the  "peaks"  and 
"valleys"  of  the  sinusoidally  varying  electric  field  and  the  observed  phenomenon  of 
the  plasma  jumping  from  peak  to  peak.  As  the  plasma  approaches  a  peak  of  the 
electromagnetic  power  the  propagation  velocity  increases.  The  plasma  is  most 
stable  at  the  peak,  accounting  for  the  slower  propagation  as  the  plasma  moves  off 
the  peak  of  the  electric  field.  The  figure  also  shows  that  the  plasma  propagation 


Figure  3.  Plasma  propagation  velocity  time  histories  for  helium  at  one 
atmosphere  and  1000  W  input  power. 
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using  the  waveguide  corrections  is  lower  than  that  of  the  planar  propagating  wave. 
Also,  the  two  cases  show  different  wavelengths  of  oscillation  and  amplitudes. 

Figure  4  shows  a  representative  electric  field  distribution  for  the  waveguide 
case  at  a  time  of  5.0  seconds  corresponding  to  a  maximum  in  the  propagation 
velocity.  As  is  evident,  the  input  field  is  totally  attenuated  in  the  results  over  a 
distance  of  approximately  10  cm.  The  attenuation  is  achieved  by  both  absorption 
and  reflection  which  will  be  discussed  subsequently.  Negligible  transmission  imder 
these  conditions  was  also  reported  in  the  experiments  of  Batenin,  et  al.  (1976b). 


Figure  4.  Representative  electric  field  distribution  for  helium  at  one 
atmosphere  and  1000  W  input  power  for  the  waveguide  case. 


Figure  5  shows  the  dependence  of  averaged  plasma  propagation  velocity  on 
input  power  for  the  two  cases  of  the  present  results  as  well  as  the  experimental  and 
theoretical  results  of  Mueller  and  Micci  (1989).  Neither  model  is  able  to  reproduce 
the  pronounced  dependence  on  input  power  observed.  The  jump  in  the 
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experimental  data  near  1500  W  may  be  attributable  to  a  different  propagation  mode 
due  to  resonant  radiation  (Batenin,  et  al.,  1973;  Mueller  and  Micci,  1989),  but  the 
lower  power  values  do  not  agree  or  show  a  similar  power  dependence.  The  one 
dimensional  model  developed  by  Mueller  and  Micd  (1988)  for  a  planar  propagating 
wave  shows  much  lower  plasma  propagation  velocities  in  comparison  to  the  present 
results,  which  vary  from  approximately  0.2  to  2.0  cm/sec  over  the  same  power 
range.  Their  results  also  show  less  dependence  on  input  power.  This  discrepancy 
may  be  partially  due  to  differing  values  for  the  properties  of  the  fluid  such  as 
conductivity  and  permittivity.  The  importance  of  the  properties  shall  be  discussed 
later.  Their  model  also  approximates  the  reflection  of  the  electric  field  from  the 
plasma  wave.  Continuous  reflection,  included  in  the  present  formulation,  allows  for 
spatially  continuous  reabsorption  of  electric  field  waves  reflected  from  within  the 
initial  regions  of  the  plasma,  yielding  steeper  temperature  gradients. 


Power  (W) 

Figure  5.  Plasma  propagation  velocity  for  helium  at  1  atmosphere 
pressure  and  500  to  3000  W  input  power. 
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Figure  6  shows  the  dependence  of  averaged  plasma  propagation  velocity  on 
the  static  pressure  within  the  chamber,  or  waveguide,  along  with  the  experimental 
results  of  Batenin,  et  al.  (1977)  and  Mueller  and  Micci  (1989).  Both  the  planar  and 
waveguide  cases  fail  to  show  the  pronounced  pressure  dependence  of  the 
experimental  data.  The  calculated  pressure  dependence  of  the  plasma  propagation 
velocity  is  approximately  p'^.  The  experimental  data,  however,  are  better 
represented  by  an  exponential  pressure  dependence,  i.e.  ,  Uprop  ~  exp(ap).  The 
exponential  factor,  a,  is  approximately  -  4.7  for  Batenin,  et  al.  and  approximately 
-  6.9  for  Mueller  and  Micci.  The  exponential  pressure  dependence  of  the 
experimental  data  further  suggests  that  the  assvimption  of  LTE  is  incorrect,  even  at 
these  relatively  high  pressures.  Non-LTE  behavior  in  relation  to  microwave 
plasmas  has  been  indicated  in  the  spectroscopic  data  of  Batenin,  et  al.  (1977, 1976c) 
and  Mueller  and  Micci  (1989).  In  this  latter  study,  electron  temperatures  were 
measured  to  be  approximately  13,000  K  compared  to  calculated  ion  and  neutral 
particle  temperatures  of  approximately  8,000  K  for  input  powers  of  500  to  3000  W. 
Recently,  non-LTE  behavior  has  been  discussed  theoretically  and  experimentally  in 
the  region  downstream  of  a  radio-frequency  discharge  in  argon  (Owano,  et  al., 
1991).  The  observed  electron  concentrations  at  one  atmosphere  and  9000  K  were 
substantially  higher  than  equilibrium  levels. 

More  insight  as  to  the  extent  of  non-equilibrium  behavior  is  obtained  from  a 
theoretical  analysis  which  follows  the  spirit  of  those  done  by  Raizer  but  with  further 
simplifications.  The  Ohmic  heating  term  is  approximated  by  the  maximum  electric 
field  amplitude  and  the  values  of  electrical  and  thermal  conductivity  are  based  on 
the  plasma  temperature.  A  simple  scaling  analysis  of  the  energy  equation  using  the 
assumption  of  total  absorption  leads  to  the  following  relation: 
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For  an  input  power  of  1 000  W  and  helium  at  one  atmosphere  and  300  K,  the  exit 
temperature  is  8000  K.  The  value  of  Emax*  which  relates  to  the  parameter  A  in  the 
theoretical  formulations,  is  approximately  1.75  x  10^  V/m  for  a  planar  propagating 
plasma  wave.  The  density  at  the  inlet  is  taken  to  be  .16  kg/m^  and  the  value  of  Cp  is 
approximately  5200  J/kg/K.  At  8000  K,  the  electrical  conductivity  is  0.65  mho/m 
and  the  thermal  conductivity  is  1.85  W/m/K.  These  values  3rield  a  propagation 
velocity  of  approximately  0.2  m/s,  which  is  in  order-of-magnitude  agreement  with 
the  present  calculations. 

Due  to  the  present  assumptions,  the  thermal  conductivity  is  independent  of 
pressure,  and  the  electrical  conductivity  has  a  weak  pressure  dependence  within  the 
range  0.4  to  1.0  atmospheres.  The  theoretical  dependence  of  the  plasma 
propagation  is,  therefore,  mainly  due  to  the  density,  giving  a  pressure  dependence  of 
approximately  p'^  for  equilibrixim  conditions.  The  experimental  data  suggests  a 
higher  pressure  dependence  for  the  electrical  conductivity.  Furthermore,  it  can  be 
inferred  that  if  the  thermal  propagation  mechanism  is  valid,  then  the  electrical 
conductivity,  and  thus  the  electron  number  density,  is  an  order  of  magnitude  higher 
than  the  LTE  value. 

Although  a  non-LTE  model  should  be  incorporated  first,  another  possibility  of 
the  anomalous  behavior  is  that  the  thermal  propagation  mechanism  may  not  be 
valid  for  these  pressures.  Batenin,  et  al.  (1976b)  suggests  that  the  propagation 
mechanism  is  not  clearly  defined  between  heat  conduction  and  resonance-radiation 
diffusion.  If  radiation  were  important,  the  thermal  conductivity  would  increase, 
possibly  giving  better  quantitative  results. 
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Figure  6.  Plasma  propagation  velocity  for  heUum  at  1081.7  W  and  0.5  to 
1.0  atmosphere  chamber  pressure. 

Figure  7  shows  the  averaged  chamber,  or  waveguide,  exit  velocity  over  a 
range  of  input  power  levels  for  the  two  cases  of  the  present  results  along  with  the 
theoretical  results  of  Mueller  and  Micci  (1989).  In  the  theoretical  work  of  Mueller 
Micci  (1989)  the  exit  velocities  were  calculated  for  both  one  and  ten 
atmospheres  of  chamber  pressure.  Qualitative  agreement  exists  between  the 
models  for  the  overall  decrease  of  exit  velocities  for  increasing  pressure.  The 
present  results,  however,  show  higher  exit  velocities  and  greater  dependence  on  the 
input  power.  Figure  8  shows  the  dependence  of  the  averaged  chamber  exit  velocity 
on  pressure  for  the  two  cases  of  the  present  results.  Note  that  the  results  for  the 
waveguide  corrected  solution  in  both  Figure  7  and  8  are  less  than  that  of  the  planar 
results,  although  each  shows  the  same  power  and  pressure  dependence. 
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Figure  7.  Chamber  exit  velocity  for  helium  at  .5  and  1  atmosphere 
pressure  and  500  to  3000  W  input  power. 


Figure  8.  Chamber  exit  velocity  for  helium  at  .5  to  1  atmosphere 
chamber  pressure. 
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Figure  9  shows  the  temperature  proiiles  of  the  propulsion  chamber 
(waveguide)  over  a  range  of  input  power  levels  for  the  waveguide  corrected  case.  All 
profiles  were  taken  at  the  same  time  level  (5  seconds).  Initial  startup  acoustic 
disturbances  were  small  and  quickly  died  away  (less  than  2  seconds  in  real  time  for 
the  lowest  power  level).  The  calculations  show  a  steep  wavefront.  For  the  500  W 
case,  the  wavefront  is  approximately  twenty  centimeters  across  for  a  temperature 
change  of  7500  K  The  planar  propagating  results  are  very  similar  to  the  case 
shown,  having  slightly  wider  plasma  wavefront  thicknesses  and  slightly  lower  peak 
temperatures.  The  values  of  plasma  wavefront  thickness  are  much  smaller  than 
predicted  by  the  model  of  Mueller  and  Micci  (1988),  especially  for  the  500  W  case. 
The  peak  temperatures  of  the  current  work  are  also  approximately  1 000  K  higher, 
again  suggesting  partial  absorption  of  the  reflected  field  as  a  source  of  the 
difference. 


Figure  9.  Temperature  profiles  for  helium  at  1  atmosphere  and  500  to 
3000  W  input  power. 
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The  absorbed  and  reflected  power  levels  in  the  present  calculations  were 
obtained  by  two  methods.  For  the  reflected  power,  the  value  of  E  x  H  was  calculated 
at  the  inlet  using  Equation  (27)  and  compared  with  the  input  value  from  Equation 
(14).  For  the  absorbed  power,  the  integral  of  the  Ohmic  heating  term,  aE2,  was 
calculated  using  Equation  (16).  The  two  values  agreed  within  three  percent  of  each 
other.  Figure  10  shows  the  time  history  of  the  percentage  of  input  power  absorbed 
by  the  plasma.  As  in  the  case  of  the  time  history  of  propagation  velocity,  the 
absorbed  power  cycles  sinusoidally  as  the  plasma  moves  through  the  regions  of  the 
peak  electric  field.  The  percent  of  power  absorbed  for  the  waveguide  case  is  less 
than  that  of  the  planar  propagating  case.  This  is  due  to  the  assumption  of  the 
sinusoidal  cross  sectional  shape  as  opposed  to  a  plane  wave  at  peak  electric  field 
strength.  Also,  note  the  longer  wavelength  and  smaller  amplitude  of  oscillation  in 
the  waveguide  case  just  as  in  the  propagation  time  history  results. 


Figure  10.  Absorbed  power  history  for  helium  at  1  atmosphere  and 
1000  W  input  power. 
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^gures  11  and  12  show  the  percentage  of  input  power  absorbed  over  a  range 
of  input  powers  and  pressures.  In  the  experimental  results  of  Mueller  and  Micci 
(1989)  for  heliiun  at  one  atmosphere,  the  power  absorption  was  within  the  middle 
60  to  lower  70  percent  range.  This  would  suggest  better  agreement  of  the  planar 
results  to  the  experimental  work  than  that  of  the  waveguide  corrected  case.  For 
helium  at  1081.7  W,  the  absorption  decreased  from  approximately  65%  at  .25 
atmospheres  to  60%  at  one  atmosphere.  Thus,  both  calculations  and  experiment 
yield  little  dependence  on  pressure,  and  overall  agreement  in  absorbed  and  reflected 
power  levels  is  considered  good. 
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Figure  11.  Absorbed  power  for  helium  at  .5  and  1  atmosphere,  500  to 
3000  W  input  power. 
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Figure  12.  Absorbed  power  for  helium  at  .5  to  1  atmosphere  chamber 
pressure. 

Further  calcxilations  have  been  done  to  ascertain  the  affects  of  a  low  speed, 
stabilizing  inlet  mass  flux.  Using  a  previous  calculation  of  the  propagation  speed  to 
determine  the  proper  mass  flux  necessary,  it  was  possible  to  stabilize  the  plasma 
near  a  location  of  maximum  power  absorption.  This  stabilization,  however,  is 
greatly  dependent  on  the  pressure  and  temperature  remaining  at  fixed  levels  and  a 
non-varying  flow  speed,  making  simple  mass  flux  control  unfeasible  for  plasma 
stabilization.  Besides  the  use  of  a  bluff  body  to  stabilize  the  plasma  in  the 
recirculating  regions  (Mueller  and  Micci,  1990),  it  is  believed  that  the  use  of  a 
simple  divergently-shaped  chamber  liner  (made  of  dielectric  material)  could  also 
stabilize  the  discharge  over  a  range  of  operating  conditions.  Although  the 
properties  of  the  flow  would  vary  somewhat,  the  plasma  would  simply  move  up  the 
diverging  duct  toward  increasing  flow  velocity,  or  down  the  duct  toward  slower 
velocities,  to  that  position  where  the  inlet  mass  flux  stabilizes  the  plasma. 
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A  computational  model  of  a  microwave  induced  plasma  proptdsion  system  is 
developed  in  one  dimension  for  the  TEio  mode  of  operation  in  a  rectangular 
waveguide.  The  propellant  used,  helium,  is  assumed  to  behave  as  a  perfect  gas  in 
local  thermodynamic  equilibriiun.  The  effects  of  radiation  are  neglected  due  to  the 
relatively  low  temperatures  and  ionization  fraction  of  the  plasma.  Also,  the  electric 
field  is  assumed  to  reach  a  steady  solution  at  each  time  interval  based  on  the  known 
frequency  of  oscillation.  Maxwell's  equations  are  solved  implicitly  using  central 
differences  for  the  complex  amplitudes  of  the  electric  field  at  each  grid  cell  center. 
The  electric  field  is  coupled  with  the  Navier-Stokes  equations  which  are  solved  by  a 
time-accurate  finite  volume  method.  This  model  sets  no  restrictions  on  the  location 
of  the  plasma.  The  time  step  is  limited  only  by  accuracy  in  resolving  features  such 
as  wave  thickness  and  propagation  speed.  The  grid  is  restricted  only  in  that  the 
cells  must  be  small  enough  to  accurately  resolve  the  plasma  wavefront.  Two 
approaches  are  developed  in  the  modeling  of  the  electric  field.  In  the  planar  case, 
the  electric  field  is  assumed  to  be  a  planar  propagating  wave.  In  the  waveguide 
case,  the  electric  field  is  assiuned  to  have  a  sinusoidal  cross  sectional  shape  and 
appropriate  corrections  on  the  phase  constant  due  to  waveguide  propagation  are 
employed. 

The  model  was  capable  of  stabilizing  the  plasma  at  a  fixed  Eudal  position  by 
adjusting  the  inlet  mass  flux.  However,  owing  to  the  eigenvalue  behavior  of  the 
mass  flux,  this  approach  to  stabilization  would  not  be  adequate  in  practice.  It  is 
suggested  that  a  simple  diverging  duct  could  be  used  for  stabilization  which  would 
hold  the  plasma  at  the  proper  location  for  countering  the  mass  flux,  even  if  the  mass 
flux  varies  to  some  degree. 
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A  simple  scaling  analysis  of  the  dependence  of  the  plasma  propagation 
velocity  on  pressure  is  calculated.  The  theoretical  analysis,  based  on  the  property 
formulations  shown  in  the  Appendix,  suggests  a  dependence  of  approximately  p*^  on 
the  propagation  velocity.  The  experimental  data,  however,  show  an  exponential 
dependence  on  pressure.  This  suggests  that  the  LTE  assumption  is  incorrect  at 
atmospheric  pressure  and  as  such,  a  non-equilibrium  thermodynamic  analysis  is 
required.  Also,  it  would  be  of  great  interest  to  have  experimental  data  at  higher 
pressures  (on  the  order  of  10  atmospheres)  to  examine  under  what  conditions,  if 
any,  the  LTE  assumption  is  valid. 

The  model  does  show  the  observed  phenomenon  of  the  jumping  of  the  plasma 
towards  the  microwave  source  due  to  the  peaks  of  electric  field  strength.  The 
percentage  of  input  power  absorbed  calculated  by  the  model  agrees  fairly  well  with 
experimental  results.  The  planar  case  shows  better  quantitative  agreement  than 
the  waveguide  case.  Neither  case  agrees  with  the  experimentally  shown  pressure 
dependence. 

Based  on  the  successful  features  of  the  one-dimensional  model,  a  two 
dimensional  formiilation  is  also  posed  for  the  axisymmetric  solution  of  a  cylindrical 
waveguide.  Although  no  results  from  the  two  dimensional  formulation  are 
presently  available,  the  model  could  be  used  to  analyze  plasmas  in  configurations 
more  representative  of  actual  propulsion  chambers.  This  would  enable  an 
examination  of  the  thermal  effects  across  the  cross  section  and  wall  convective  heat 
transfer. 
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A.1  NOMENCLATURE 


PROPERTIES 


Property  Table  Variables: 

C+  ion  concentration 

fHe  electron  partition  function  for  the  helium  atom 

fHe"*^  electron  partition  fimction  for  the  helium  ion 

h  Planck's  constant 

Ii  first  ionization  energy 

K  Boltzman's  constant 

m«  electron  mass 

n  number  density  of  particles 

n«  number  density  of  electrons 

p  pressure 

Q  collision  cross  section 

q.  electron  charge 

T  temperature 

vth  thermal  velocity  of  particle  collisions 

e  permittivity 

Eo  permittivity  of  free  space 

V  collision  frequency 

a  electrical  conductivity 

0)  angular  microwave  frequency 

0)p  angudar  plasma  frequency 
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A.2  FORMULATION 


The  computation  of  the  electromagnetic  and  fluid  equation  systems  requires 
various  property  values  of  the  fluid  (gas)  used.  These  values  are  found  through  a 
table  look-up  using  a  simple  linear  interpolation  scheme.  The  table  was  developed 
such  that  all  properties  are  indexed  by  pressure  and  temperature. 

The  values  of  Cp,  Cv,  7,  are  found  \ising  the  Saha  relations.  As  such,  the  gas  is 
assumed  to  be  an  equilibrium  mixture  of  ideal  gas  species.  For  example,  heliiun  gas 
would  be  considered  a  mixture  of  He,  He+,  and  electrons,  all  at  the  same 
temperature.  Let  C+  represent  the  concentration  of  helium  ions  in  the  mixture. 
Then  using  the  Saha  equation,  the  concentration  of  ions  is  given  by 


(1-CJ(1+CJ 


(Al) 


where  all  units  are  in  mks.  The  electron  partition  flmctions  for  the  ion  and  neutral 
atom  have  the  basic  form 


f 


=  go  +  gi  exp 


h’cwil 
1  KT  / 


(A2) 


where  the  values  for  go,  gi,  and  wi  (cm  *i)  were  taken  from  Atomic  Energy  Levels  by 
C.E.  Moore  (1949).  The  properties,  c  p,  c  y,  Y»  and  the  density,  p,  are  then  calculated 
using  perfect  gas  relations  (Matsuzaki,  1982). 

The  values  for  viscosity  and  thermal  conductivity  were  taken  from  Tables  on 
the  Thermophysical  Properties  of  Liquids  and  Gases  (Vargaflik.  1975)  and  the  CRC 
Handbook  of  Chemistry  and  Physics  (Weast,  1986)  respectively.  This  data  does  not 
include  the  entire  temperature  range  needed,  so  a  linear  extrapolation  of  the  data 
was  done.  Also,  it  was  assumed  that  there  is  no  pressure  dependence  for  the  range 
used. 
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Since  the  permittivity  and  electrical  conductivity  of  the  gas  depend  on  the 
microwave  frequency,  the  values  of  the  plasma  frequency  squared,  0)p2,  and  the 
effective  collision  frequency,  v,  are  stored  in  the  table.  The  values  of  the 
permittivity  and  conductivity  are  then  foimd  within  the  program  by  the  formulae 

eovci^2 

0)2  +  v^  (A3) 


e  =  eo  1  - 


where  to  is  the  angular  frequency  of  the  input  microwave.  The  plasma  frequency 


squared  is  given  by 


(Op2=Beael 

Eonie 


where  ng  is  the  number  density  of  electrons,  qg  is  the  charge  of  an  electron,  e©  is  the 
permittivity  of  free  space,  and  mg  is  the  mass  of  an  electron.  The  effective  collision 
frequency  is  given  by 

v  =  v-  +  v-n  (AR) 


where  Vgj  represents  the  electron-ion  collision  frequency  and  Vgn  represents  the 
electron-neutral  atom  collision  frequency.  In  each  case,  the  collision  frequency  may 
be  expressed  as 

V  =  n  Q  v*  (A7) 


where  Q  is  the  cross  sectional  area  for  the  collision  and  vth  is  the  thermal  velocity 
which  is  assumed  to  be  a  Boltzman  distribution  and  is  thus  given  by 
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(A8) 


The  cross  section  for  electron-ion  collisions,  the  coulomb  cross  section,  can  be 
accurately  calculated  theoretically  (Jahn,  1968). 


a  = _ isl _ 

(47ceof(|K'l)" 


In 


8neo(|KT) 

V  n^Qe^ 


(A9) 


The  cross  section  for  electron-neutral  collisions,  however,  is  based  on  an  empirically 
derived  formula  taken  from  the  work  by  D.  E.  Golden  (1966),  The  value  used  in  the 
program  is  twice  that  of  the  cross  section  found  in  the  referenced  work.  This  was 
done  to  achieve  closer  pressure  dependence  of  the  properties  as  shown  in  the  works 
of  Batenin,  et  al.  (1977)  and  Mueller  and  Micci  (1989).  Following  are  some  plots  of 
the  resulting  property  values  derived  from  the  tables. 


Figure  A 1 .  Specific  heat  values  for  helium  at  .5  and  1  atmosphere 

pressure. 
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Viscosity  (x  10  *  kg/m/s) 


Figure  A2.  Thermal  conductivity  of  helium. 


Figure  A3.  Viscosity  values  for  helium. 


47 
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Figure  A4.  Permittivity  values  for  helium  at  .5  and  1  atmosphere 

pressure 


Figure  A5.  Electrical  conductivity  values  for  helium  at  .5  and  1 
atmosphere  pressure. 
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